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One of the great challenges of electronic structure theory is the quest for the exact functional of 
density functional theory. Its existence is proven, but it is a complicated multivariable functional 
that is almost impossible to conceptualize. In this paper, the asymmetric two-site Hubbard model 
is studied, which has a two-dimensional universe of density matrices. The exact functional becomes 
a simple function of two variables whose three dimensional energy landscape can be visualized 
and explored. A walk on this unique landscape, tilted to an angle dehned by the one-electron 
Hamiltonian, gives a valley whose minimum is the exact total energy. This is contrasted with 
the landscape of some approximate functionals, explaining their failure for electron transfer in the 
strongly correlated limit. We show concrete examples of pure-state density matrices that are not 
w-representable due to the underlying non-convex nature of the energy landscape. For the Hrst 
time, the exact functional is calculated for all numbers of electrons, including fractional, allowing 
the derivative discontinuity to be visualized and understood. The fundamental gap for all possible 
systems is obtained solely from the derivatives of the exact functional. 


In 1964 Hohenberg and Kohn [1] established density 
functional theory (DFT) showing that the electron den¬ 
sity, p, is all that is is necessary to determine the ex¬ 
act energy of many electron systems. However, all the 
challenge of electronic structure is then moved into an 
unknown universal functional of the density, fF{p\. For 
a wavefunction, ik,,, that is the ground-state solution of 
the Schrodinger equation with potential v, 

E, = = (^„|r + + Tr(p„u) (1) 

simply subtracting off the one-electron term, gives the 
exact Hohenberg-Kohn functional for py (ik.;, —> p„) 

E^^[Py]=Ey-Tr{pyV) = y\T + Vesl'S y) . (2) 

This procedure can be carried out for many different v, 
to obtain many points of the exact functional E^^\py\. 
A question arises of whether all possible densities are 
achievable. This is the problem of u-representability, that 
is addressed by the constrained search by Levy and Lieb 
[2, 3] following earlier work by Percus [4] 

= min(4'|r + T4e|4'). (3) 

This functional is defined for all possible densities com¬ 
ing from a A^-electron wavefunction, including those that 
are not obtainable as the ground-state solution of a 
Schrodinger equation (not u-representable). Once the 
exact functional is known, the total energy is obtained 
by minimization only over densities, 

Ey[p] = min {F[p\ + Tr(pu)} . (4) 

p 

The exact functional of the first-order density matrix, 7, 
can be derived [2, u] 

i^Levy[^] = min(vl/|p,,|vl/), (5) 

W—>-7 


and used similarly, where the kinetic energy term is now 
a known linear functional of 7 

Ey[p\ = min{^[7] -b Tr(T7) -f Tr(u7)} . (6) 

7 

In this Letter, the nature of the exact first-order den¬ 
sity matrix functional is revealed by considering the 
asymmetric two-site Hubbard model. In this universe, 
the fundamental equations are tractable and the exact 
functional becomes a visualizable three dimensional en¬ 
ergy landscape in the space of density matrices. We 
demonstrate how this one universal landscape gives the 
exact energy of all possible systems, for all numbers of 
electrons including fractional. This connected view of 
the functional for all density matrices makes clear the 
reasons for the failure of approximate functionals, and 
allows us to answer the questions of whether there are 
density-matrices which are not u-representable and also 
how the derivatives of the exact functional give the fun¬ 
damental gap. 

The asymmetric two-site Hubbard [6] model describes 
interacting electrons on a lattice of two sites that con¬ 
tains the physics of electron transfer and has even re¬ 
cently been experimentally described using two ultracold 
fermionic atoms [7]. It has the Hamiltonian 

H — t ^ 4“ ^2cT^tcr^ 4" U ^ T ^ 

(7 i i(T 

( 7 ) 

where the site index i = 1,2, spin index a = a, jS and the 
number operator is hi^ = c|g.Cio-- There has been recent 
work on the exact functional in this model from Fuks et al 
[8, 9], Carrascal et al [10], Pastor and coworkers [11, 12[, 
Requist et al [13[, and in other systems [14-16[. 

The parameters that define a particular model are the 
hopping between the sites, t, on-site energies ei/e2 and 
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from FCI caclulations 

Figure 1: Energy landscape of the exact functional (a) ^’'^'*''^[ 7 ] for all allowable density matrices of the two site Hubbard model, 
(b) The one electron term, y.u, for t = 0.1 and Ae = 0.9, which is purely a flat plane, (c) Illustration of the minimization of 
the exact functional adding on the same 'y.v term to give the FCI energy and density matrix, |} . (d) ^’'^*''^[ 7 ] and 
6552 points of {^’^^[ 7 ^'^^], 7 ^^^^} that show the subtracting the one electron term (Eq. 11) at 7 ^^^^ for many different v. 


the electron-electron repulsion penalty due to double oc¬ 
cupation of a site, U. The physics is completely deter¬ 
mined by Ae = 61 — 62 and the ratio U/t, therefore, in 
this work U is fixed at 1 and t and Ae are the chosen 
variables. The kinetic and on-site potential part of the 
Hamiltonian, which together we denote as v, is a real 
symmetric 2 x 2 matrix defined by parameters t and Ae 


V = 


Ae/2 -t \ 
-t -Ae/2 ) 


( 8 ) 


and the 2 x 2 density matrix, 'jij = is 


from an exact diagonalization full configuration interac¬ 
tion (FCI) calculation with one-electron Hamiltonian v, 
the Hohenberg-Kohn functional is given by 

=£;FCi^27i2f-7iiAe/2 + (2-7n)Ae/2. (11) 


The second way is the constrained search over real singlet 
wavefunctions 




^ [^((/icr(/)2/3) + A{(j)2a(j)iP)] 

+bA{(j)ia:(j)iP) + cA{(j)2a(j)2l3) 


( 12 ) 


7 = 


711 712 \ 

712 (2 - 7ii) / 


(9) 


leading to a total energy for real density matrices 


Ey = - 27 i 2 t-b 7 iiAe /2 - (2 - 7 ii)Ae/ 2 -I-F’[ 7 ]( 10 ) 


The exact functional can be obtained and understood 
from different perspectives. Firstly, for any 7 that comes 


which can be simplified to an expression (see Refs. [10, 
12] and supplementary information (SI) for more details) 

7?2 (1 - \/l - 7?2 - [711 - 1 ]^) + 2[7ii - 1 ]^ 
2 ( 7?2 + [ 7 ii - 1 ]^) 

(13) 

Thirdly, it can be viewed as the exact functional in den¬ 
sity matrix functional theory for two electrons. From the 
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Figure 2: Entire landscape of ^[ 7 ] for the exact functional, and three approximate density-matrix functionals (see supplementary 
information for more details). The minimizing values {T’[ 7 t)], 7 t)} for three lines of v (—2 < Ae < 2) with t = 1,0.2,0.05 are 
plotted on the surfaces in purple, green and red, respectively. The central plots show the failure of approximate functionals to 
correctly describe the electron transfer (711 vs Ae) as t approaches the strongly correlated limit (see supplementary animations). 


work of Lowdin and Shull in 1956 [17] using the natural 
orbitals |a) and \b) {\p) = X]i=i 2 their 

occupation numbers Ua and Ub that diagonalize 7 , it can 
be derived that 

^^^[ 7 ] = ^na{aa\aa) + ^nb{bb\bb) - ^yn^{aa\bb) (14) 

where the two-electron integral is {pp\qq) = 
UYl,i=i 2 ^pi^qi- This gives exact agreement with 
the constrained search expression, Eq. (13) and has been 
utilized in functionals such as the AGP natural orbital 
functional [18, 19[ and PNOF5 [20[ (see SI). There 
are two further possible routes to the exact functional 
(details in the SI): the extension over pure-state wave- 
functions to complex, and the Lieb maximization[3[, 

^Levyj-^j jg gjjown in Fig. la. for the allowable density 
matrices (711 — 1) -I- 7^2 <1- It is represented as a 
unique surface of hills and a valley in a bowl type shape, 
with a channel through the centre (at 711 = 1 ) and hills 
on both sides (reaching 1 at 712 = 0). This defines the 
energy landscape that maps every possible system to its 
corresponding exact energy. 

The exact functional is an energy landscape with only 


one minimum, so how does it give rise to all possible FCI 
energies? This can be pictured in a very physical manner 
by considering a walk on this landscape, placed upon a 
flat surface tilted to the angle given by the one-electron 
potential, which gives a valley whose minimum equals 
exactly the FCI solution. Fig. lb shows the one electron 
term for a particular v, defined by t = 0.1 and Ae = 0.9, 
and Fig. Ic shows the addition of this with the exact 
functional, F'^®''^[ 7 ] -I- 7 .U, whose minimum is at the FCI 
energy, and FCI density matrix, 7 „. This holds 

for every possible v. Thus, once the exact functional is 
known, it gives the exact solution of any system by means 
of an almost trivial calculation. 

We have performed a large number of FCI calcula¬ 
tions varying the two free parameters, —10 < t < 10 
and —10 < Ae < 10. Fig. Id illustrates the result of 
over 6000 FCI calculations subtracting off the one elec¬ 
tron term, 7 .U, to give the ^^^[ 7 ] of Eq. (11). Every 
single light blue dot, representing many F^^[ 7 j,], lies on 
the surface of F^®''y[ 7 ]. However, the one-particle den¬ 
sity matrices 7 „ that result from all these FCI calcula¬ 
tions cover only a small fraction of the space (seen as the 
black dots projected onto the base of the plot with more 
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F'“®''''['y] for Number of electrons = 1.5 
Line of symmetric Hubbard models for N=1.5 
=0.5 0.5E,^C'[2]-y.hi 



Figure 3: The exact functionals of Eqs (17) and (18) for 
= 1.5 electrons, which gives back the exact energy of every 
system with 1.5 electrons. 


Fn+sM 

Line of symmetric Hubbard models for N=1.5 



Figure 4: The exact functionals of Eqs. (17) and (18) for all 
numbers of electrons (1 < < 3) in the symmetric two site 

Hubbard model. 


details in SI). The rest of the density matrices are not 
w-representable, even though they are Al-representable. 
From the perspective of the exact functional, it is clear 
why these density matrices can never be found, as they 
correspond to the hills of the surface where the [7] 

lies inside a convex containing surface (see SI). Addition 
of the one electron interaction term, which is purely lin¬ 
ear in the variables 711 and 712, as pictured in Fig. lb, 
means that these points can never be minima, and hence 
cannot be a FCI solution. In terms of the functional it 
corresponds to where the second derivatives of the func¬ 
tional are no longer positive definite as seen by the a 
negative lowest eigenvalue of the Hessian matrix of sec¬ 
ond derivatives, Hij = 9777^77, (see SI). It should also be 
noted that the lowest energy wavefunctions of the non- 
u-representable density matrices cannot be written in a 
Gutzwiller form [21] (see SI). The non-u-representable re¬ 
gion highlights the key distinction between ^'^^''^[7] de¬ 
rived from pure-state wavefunctions, which can be con¬ 
cave, versus the F^'®’^[7] functional derived from ensem¬ 
bles by a Legendre-Fenchel transform, which is proven to 
be everywhere convex [3]. 

The derivatives of the functional (expressions in SI) 
satisfy the Euler equation and give the one-electron 
Hamiltonian needed, 

„ =-v + C. (15) 

07 

Now, consider the physics of electron transfer, by vary¬ 
ing Ae, from the weakly correlated {U/t = 1) to strongly 
correlated {U/t = 20) regimes as depicted in Fig. 2. 
Correctly describing this electron transfer in the strongly 
correlated regimes is one of the great challenges of elec¬ 
tronic structure, as demonstrated in Fig. 2 by the failure 
of approximate density matrix functionals such as Muller 
[22] and Power functionals [23] . The approximate func¬ 
tionals do not correctly describe the entire landscape and 



thus completely fail to describe electron transfer (see an¬ 
imation in SI). This is related to the complete failure 
of all currently used density functionals for the electron 
transfer in a two-electron molecular type challenge (see 
HZ{2HofRef. [24]). 

The exact functional can be calculated for all numbers 
of electrons {0 < N < 4); the integer parts are trivial 
and given in the SI. For non-integer numbers of elec¬ 
trons, the functional is constructed using the Perdew, 
Parr, Levy and Balduz (PPLB)[25[ ensemble extension 
to search over many-electron density matrices 

Piv+i = co|4'o)(4'o|-F ci|4'i)(4'i|-F C2|4'2)(4'2| 

Tesla's) («'3| +C4|«'4)(^'4| (16) 

with Ci = 1 and Ci.i = N + S (0 < 5 < 1). Thus, 
we explicitly construct the fractional extension 

TV+afy] = min Tr^Ar+^Pee], (17) 

rjv+5->'7 

where, unlike PPLB, we have not assumed convexity 
of the energy versus N. That is, rather than using 
PAT+a = CAr|'I'Ar)(4'Ar| -F cat-i-i|' pAf+i)('I'Ar+i|, we explic¬ 
itly search over ensembles of all A-electron wavefunctions 
{N = 0,1,2, 3, and 4) as in Eq. (16) (see SI). 

Fig. 3 shows the extension of the exact functional to 
fractional numbers of electrons for N + S = 1.5. We 
obtained FAr+afy] for all the possible density matrices, 
where the minimum is actually given only by the combi¬ 
nation of N and A -F 1 (see supplementary information 
for more details). We also find that all the appropriate 
ensembles of FCI energies subtracting off the one electron 
term using the ensemble of density matrices, 

fHK,[u] = {1-S)E^^\N]+6E^^\N + 1] 

- [(l-<^)7f+<57^“^^](18) 

lie perfectly on the functional surface for all values of v 
and 6. Additionally, just like for integer electrons, a walk 
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on this surface tilted to the angle of any one-electron 
potential (analogously to Fig. Ic) gives a minimum point 
that exactly agrees with the ensemble FCI energy. 

The knowledge of the exact functional for fractional 
numbers of electrons connects to the band-gap problem. 
This is the question of whether the fundamental gap, de¬ 
fined as the difference of the ionization energy and elec¬ 
tron affinity, can be given by the derivatives of the exact 
functional. For simplicity, consider only the symmetric 
Hubbard dimer with different numbers of electrons. In 
Fig. 4 the exact functional is shown for 1 < iV < 3, 
along with several points of the ensemble with 

V = {—1 < t < l,Ae = 0} (see also animations in SI). 
For every v, traces out a straight line versus 

particle number with a clear derivative discontinuity at 
N = 2, hence the derivatives of the exact functional give 
the contribution to the fundamental gap 


dF[^] 


dN. 


dF[j] 


dN_ 


= F[7^+1] + ^[7^"-'] - 2F[j‘^]. 


iV-ll 


..iVl 


(19) 

If there is no discontinuity in the density matrix, which 
is the case of a Mott insulator, the entirety of the fun¬ 
damental gap is given by the exact functional (Eq. 19). 
This is illustrated as the green line in Fig. 4 for the 
symmetric Hubbard model with t = 0 and 1 < iV < 3, 
and has a direct correspondence to the gap of infinitely 
stretched H2 [26] . Nevertheless, most systems have a dis¬ 
continuity in the density matrix, ^ 

giving rise to a discontinuous derivative even for the one 
electron term, which is an entirely smooth flat plane. 
However, the direction in which 7 changes upon electron 
addition or removal is already determined by derivatives 
of F whilst keeping the derivative in the direction of fixed 
N to be constant 


7Af±i = 7./V + 


^7 


SN± 


dF I 
O'Y \ N 


( 20 ) 


Hence, the fundamental gap is solely determined by the 
derivatives of the functional itself. 


Gap[7N] = 


dF[j] 


dj 


dFb] 


dN+ sfit] I 

St- In 

9A^-i_ ^ 1 

d'y \ N 

N 

I c)F['y] 


dj 


5E[7] 

1 dN_ 

8F[7] I 

9~l In 

dN_ 

dF 1 

liv 

^7 


N , 


Overall, it is amazing to have a universe that turns any 
question about the exact functional into simple move¬ 
ments of a three-dimensional energy landscape. Walks on 
this landscape and its valley and hills correspond to im¬ 
portant physical concepts such as the exact energies of ev¬ 
ery possible system and domains of non-u-representable 
density matrices. Furthermore, in the direction of chang¬ 
ing particle number there is a continuous surface that has 


a derivative discontinuity at the integers, giving all pos¬ 
sible fundamental gaps, including Mott insulators. The 
whole landscape of the exact functional is itself an infinite 
number of exact constraints, such that any approxima¬ 
tion must approach and be mathematically proximal to 
it for the entire universe. It is this connected view of the 
exact functional for a family of densities in a global land¬ 
scape that truly highlights a path for the improvement 
of approximate functionals. 
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A. Derivation of Eq. (13) 

To derive the exact functional, 

F^^^^ = min(^'|V;e|«') 

'L—>-7 


(SI) 


consider the minimization over real singlet wavefunctions 


+ A{(l)20i<fil3)] 

+bA{4>ia4>il3) + cA{(j)2a(j)2l8)- 


(S2) 


in terms of the parameters a, b and c along with the 
normalization = 1 and the elements of 

the density-matrix 7^- = giving 711 = 

2b^ -\- of and 712 = {ba + ac). The two-electron en¬ 
ergy comes only from the (11|11) and (22|22) integrals, 
which are U, as all other integrals are 0, so only the sec¬ 
ond determinant with itself and the third determinant 
with itself contribute, giving 


Fm = U{h'^+ c^) = U{l-a'^) 
It is also satisfied that 

711 - 1 = 5^ - c^. 
Therefore, using 712 gives 


{h + c) = 


712 

\/2a 


and combining with Eq. (S4) leads to 

(711 - 1 ) y/2a 


{h-c) = 


712 


Now, square Eqs. (S5) and (S6), to give 




and 


{b-cf = 


(711 - 1 )^ 20 ^ 


I12 


(53) 

(54) 

(55) 

(56) 

(57) 

(58) 



Figure SI: The plane of all possible density matrices illus¬ 
trating the non-w-representability of many of the allowable 7 . 
a) The second derivatives of the exact functional showing the 
points where the lowest hessian eigenvalue is < 0 from Eqs 
S14-S16 and b) the density matrices, 7 , achieved in 6552 FCI 
calculations for —10 < t < 10 and —10 < Ae < 10. 


Adding these two has the result 


2b^ + 2C = ^ , 

2a^ 

Using the normalization, gives 


7?2 , (7ii-l)^2a2 


7?2 


2^ 1I2 , (7ii-l)"2a2 


(2-2a^) = ^ 


7?2 


(S9) 


(SIO) 


(Sll) 


which leads to a quadratic equation for of 

with solution 

„ 7^2 (1 ± v'1-(7ii-1)^-7?2) 

a2 = ^ ^ - -g- - -A. (S12) 

2[(7 ii - 1)' + ih] 

Taking the plus combination gives the lowest energy 


E = l-a 
= 1 - 


I12 (1 + - (711 - 1 )^ - 7 ? 2 ) 


2[(7 ii - 1)^ +7?2] 

2[(7ii - 1)^ + 7?2] - 7?2 (1 + Vl - (711 - 1)^ -7?2 
2[(7ii-1)2+7?2] 
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2(711 - 1)^ + ill (l - \/l - (711 - 1)^ - 7^2) 


2[(7ii- 1)^+7?2] 
This agrees with Eq. (13) of the paper. 


The derivatives of the exact functional can be evalu- 
(S13)ated analytically and are used in Fig. SI. 


dE _ 4(711 - 1 ) + (711 - l)7i2/\/l - (711 - 1 )^ - 7i 2 f) (2(711 1 ) +712(1 ^/^ Ip ^)) 

^711 2 ((711 - 1)2 - 7^2) _ 1)2 + ^2^]2 

nd 

QE ^?2/\/l - (711 - 1 )^ - 1I2 + 2712 (1 - yi - (711 - 1 )" - 7 p) 2712(711 - 1 )^ + ypx/i - (711 - 1 )" - 1I2 

^712 2 (( 7 ii - 1)2 + 722) [(.^J^_l) 2 +^ 2 ^j 2 


d^E 


- 4(711 - 1 )" ((711 - 1 )" + 71=2) ( , , ^"2 + 4 ] + 8(711 - 1)2 (2(711 - 1)2 - 7?2 ( 


-7?1 + 2711 - 7=2 - 1 


)) 


((711 - 1)2 + 7 ^ 2 )^ ( 4 - 


(-^?1+2 t 11-T?; 


2 ((711 - 1)2 + 7 ^ 2 ) 

p 72 ^ - 2 ((711 - 1)2 + 7 ^ 2 ) ( 2(711 - 1)2 - 7?2 (\/- 7 ?i + 2711 - 7?2 - 1 
2 ((711 - 1)2 + 7 ^ 2 )^ 


+ 814 ) 


d^E 

97110712 


(711 - 1)712 (-27^1 + 1275^ + 272^ (i0a/- 7P + 2711 -7^2 - 97^2 + 1) - llli (2 a/-7?i + 27ii - 7?2 “ ^7^2 + ^) ) 

2 (-7^ + 2711 - 7^2) (7?! - 2711 + 7?2 + 1)^ 

(711 - 1)712 (7^2 (-27?2 (2\/-7?i +2711 -7^2 + 3 ^ + 4 a/- 7 P +2711 -7^2 + 7^2 + l) ) 

2 (- 7 ?i + 2711 - 7^2)^^^ (7?! - 2711 + 7?2 + 1)^ 

(711 - 1)712 (+ 7^1 (4 (\/- 7 ?i + 2711 - 7?2 - 6 ^ - 37 ^ 2 ^ - 473^ ( 4 a/- 7 ?i + 2711 - 7?2 “ 37?2 “ ^ 


92 E 


2 (-7?! + 2711 - 7^2)^^^ (7?! - 2711 + 7^2 + 1)" 
-47?2((7ii-1)^+7?2) ( ' / 2 T -2'-V-Pi+27ii-7?2 + 2) 

VV-+"i+ 24 ii -^?2 J 


(S 15 ) 


87? 


2 ((711 - 1)^ +7^2) 

(2(711 - 1 )^ - 7?2 (\/- 7 ?i + 2711 -7?2 - 1)) - 2 ((711 - 1 )^ +7^2) (2(711 - 1 )^ -7?2 (\/- 7 ?i +2711 -7 i2 “ l)) 

2 ((711 - 1)2 +7^2)^ 


((711 - 1 )^+7^2)^ ( / 2 ~ ^ V'-'ln + 2711 - 7 i 2 + — 2,2^"^ 2 -( 3/2 +2) 

\V“^ii+ ^““^12 (-+i+2^ii-+2) J 


2 ((711 - 1)^ +712)" 


(S 16 ) 


For the discussion of u-representability, there are two 
common counterexamples: the first is a one-electron 
density with a certain type of cusp, given by Englisch 
and Englisch [1]; the other is a spherical p density re¬ 
lated to a degeneracy that cannot be given by a single 
wavefunction['2]. The non-u-representable density matri¬ 
ces shown here are very different to these two examples 


and are only due to the nature of the energy surface of 
the exact functional as shown in Fig. S2. 

B. Derivation of Lowdin-Shull for Hubbard model 

Lowdin and Shull (LS) showed that the natural or¬ 
bitals, (fik, that diagonalize the density matrix and wave- 






















































Figure S2: Two lines of that illustrate non-w- 

representable density matrices, due to the non convexity of 
the surface along the given line. 


function for two electrons are the same 


^(r,!-') = ^Cfc<()fc(r)<()fc(r') 
k 

(S17) 

7 ( 1 ', r') = '^nk(t)k{-r)(t)k{F) 

(S18) 


k 


where = 2c^. 

2 

='^ct {2hu + {ii\ii)} + 2ciC2{ll\22) (S19) 

i=l 

For two basis functions the minimum energy wavefunc- 
tion comes from the coefficients of ci and C 2 having op¬ 
posite signs, Cl = y/ ni/2 and C 2 = —\fn 2 l 2 . Substitut¬ 
ing this into the energy expression for the wavefunction 
gives an expression in terms of the natural orbitals and 
the natural orbital occupation numbers, rifc, 

^^^[7] = ]^na{aa\aa) + ]^ni,{hb\hb) - y/n^{aa\bb). 

(S20) 


There has been some recent interest in natural orbitals 
[3] and natural orbital functionals that, for two electron 
systems, must reduce to the Lowdin-Shull expression if 
they are to be exact, for example the PNOF5 functional 
[4-6], 


The eigenvalues of the density matrix 7 = 

7 11 712 

712 (2 - 711) ' 


(711 - n)((2 - 7ii - n) - 7^2 = 0 (S21) 

- 2n + 2711 - 7 n - 71^2 = 0 (S22) 


n± = (^-2 ± ^4- 4(711 -1)2+ 4722^/2 (S23) 

n± = ria/b = 1 ± (711 - 1)^ + 7 i2 - (S24) 


The {pp\qq) integrals are in the natural orbital basis 
and the coefficients of the natural orbitals (Cpi) are 
found by substituting in the natural orbital numbers e.g. 
(711 - np) Cpi -1-712Cp2 = 0 or (also using C±i = Cf^a/bp) 


C±i 


(711 - 1) ± 1/(711 - 1)2 -h 72^ 

712 


C ±2 and Cly+Cl^ 


(S25) 

So overall, C|.i = a^±l (712 + o-±) and C '^2 = 712 /(a± + 
712 ) and hence 


= lria{Ct,+C^,2)U+lnb{Ct,+Ct2)U 

-VFJibiClyCi, + Cl2Cl2)U. (S26) 


For convenience, replace r = (711 — 1) and S = 
y/r'^ + 7 i 2 , to obtain the following expression 


^LS _ 


(l + S) 


(r -I- 5)2 


7?2 + (^ + ^)^ 


-h 


112 


-7/1-52 


(r -I- 5)2 


il2 + ir + sy 
(r-Sr 


( 1 - 5 ) 


7?2 


(r- 5)2 


7 h + ir-sy 

7?2 


7?2 


ih + ir-sy 


1I2 + {r + 5)2 7^2 -t (r - 5)2 722 + {r + 5)2 + {r - 5)2 


(S27) 


This equation could be simplified further but we have 
checked, by numerical evaluation with Fortran code, that 
it gives identical results to Eq. (13). 


C. Complex 

The constrained search 'F —?► 7 can be expanded over 
complex wavefunctions where the parameters, a, 6 , c, in 
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from complex 'F 



Figure S3: The exact functional allowing the wavefunction to 
be complex in the Levy search. 


the wavefunction 


-^compiexH lower than Eq. (13) the solutions have a 
current and this may give a connection to the exact func¬ 
tional in current DFT (CDFT) [ 8 , 9]. 


D. Lieb maximization 


f'-'®*’ 6552 FCI calculations 





^ [A{4>iOi4>2l3) + A{4>20i4>ili)] 
+bA{(j)ia(j)ij3) + cA{(j)2Q-(j)2li) 


(S28) 


are allowed to be complex 


a = ar + iai 
b = br + ibi 
C = Cj’ “1“ %Ci 

In terms of these parameters there are the following con¬ 
straints: 


1 = al + aj + bl + bj + cl + cj 

Til = 2a^ -|- 2a^ b^ b^ 

^( 712 ) = V2{arbr + Qibi -I- brCr + hd) 

The imaginary part 3 ( 712 ) can be anything as it does 
not enter the energy expression. A fourth constraint can 
be included if the overall phase of the wavefunction is set 
to zero. 

We now carry out a search over all possible wavefunc- 
tions minimizing E and a given 711 and 5 i( 7 i 2 ), which 
gives Fig. S3. We do this by an explicit grid search 
over the two remaining variables for each 711,712 that 
is specified. The resulting energy functional gives the 
same result as the Hubbard expression Eq. (13) for 
all density matrices except the non-u-representable set. 
For all possible FCI density matrices it is, of course, in 
agreement with F^^[ 7 „]. For the non-w-representable 
set, ^(.ompiexlT] lower in energy, though this does 

not change any physics as these points can never be 
minima of any Hamiltonian. In this case, the func¬ 
tional numerically agrees with the ensemble functional 
considered by Saubenere and Pastor[7] given by a den¬ 
sity matrix that is an ensemble of two wavefunctions 
r = a|T'a)(4'al(T'bj. It should be noted that when 


Figure S4: Functional from Lieb maximization using 

6552 FCI calculations 

Another way to to calculate a bound for the functional 
is to perform the Lieb maximization[10], 

F^^^'^[-f]=sup{E,-7.v} (S29) 

V 

which is a supremum (a smallest upper bound which for 
any finite set would just be a maximum) on the set of v. 
This means for a finite set it would actually be a lower 
bound to the true minimum ^'^“*’[ 7 ] < ^'^^''^[ 7 ]. The 
Lieb maximization is carried out using 6552 FCI calcu¬ 
lations for V, with —10 < t < 10 and —10 < Ae < 10. 
Over a grid of density matrices, we compare directly with 
-^complex from complcx wavcfunctions as in the region of 
non-w-representable densities it is closest to the complex 
or ensemble form. Carrying out the maximization of Eq. 
(S29) gives the results in the left hand side of Fig. S4 and 
the difference to is shown in the right-hand side. 

This difference is small and negative which illustrates 
that the Lieb maximization only gives a lower bound to 
the true functional that in this case is known exactly. 
Obviously, with more and more FCI calculations 
would approach closer to the correct result. The ^^“*’[ 7 ] 
should not be used in minimizations in the same way 
as as it is a lower bound rather than an upper 

bound. Finally it should be noted that ^'"“'^[ 7 ] is every¬ 
where convex by construction and cannot, for example, 
contribute to the discussion on u-representability. 


E. Approximate Density Matrix Functionals 

We consider various approximate density matrix 
functionals including Hartree-Fock as a density ma¬ 
trix functional, Muller[ll], Power [12]. Here the 
value of the natural orbital occupation numbers 0 < 
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rii < 2 and the two-electron integrals {pq\rs) = 

f f (/>*(r)(/>r(r)Vee(r,r')(/>*(r')(/>s(r')drdr' which in the 
asymmetric two-site Hubbard model just work out to be 
(pqjrs) = X]i=i 2 ^pi^qi^riCsi in terms of the orbitals 
coefficients Cpi (|p) = C'p*4|vac)) 


F 


Hartree—Fock 


= \nrnj{ij\ij) - ^ninj{ii\jj) 


If we consider all possible values of ci and — 1 < g < 1 
we get the following density matrices and 


F[dJ 


GWAi 


(^GWA|y^^|^GWA^ 
^VJ/GWA|vJ/GWA^ ■ 


For other values of jgl > 1 the wavefunction is no longer 
a ground state wavefunction. 


^Muller ^ ^mrijiijlij) - ^^n,nj{ii\jj) 


pPowev ^ ^mnj{ij\ij) - ^{ninj)°‘{ii\jj) 

In the paper we use a value a = 0.675 that has recently 
been used for Mott insulators 


F. Gutzwiller approximate wavefunction 

The Gutzwiller wavefunction [13] is a parametrized 
wavefunction of the form 

^ [^((/)iQ;(/)2/3)+-4.(^2Q;(/>1;8)] 

+g [A{(j)ia(j)iP) + ^((/)2Q!(/)2/3)] 

When 5 = 1 it is the Hartree-Fock wavefunction for or¬ 
bitals (j) = + <^2)- The basic idea is that in an H2 

like system as 0 it goes to the Heitler-London wave- 
function. In the asymmetric two-site Hubbard model we 
consider an orbital of the form 4> = ci(j)i + -^/l — c\(j)2 
and a Gutzwiller wavefunction 

Gutzwiller 0<g<1 + 





G. Functional for = 0,1, 2, 3 and 4 

The functional is calculated for different integer num¬ 
bers of electrons {N = 0,1,2,3 and 4), where the trace 
of the density matrix yn -|- 722 = N. At iV = 0, 
F[7] = 0 and there is only one allowed density matrix 
7ii = 712 = 0. For = 1, F[7] = 0 as there is no 
electron-electron interaction, however, the allowable den¬ 
sity matrices from a pure state wavefunction are now de¬ 
fined by a circle 712 = (711 — 0.5)^ — 0.5^. Inside this 

circle are ensemble-Wrepresentable density matrices but 
they cannot come from a pure-state wavefunction. For 
N = 2, F'[7] is that of Eq. (13). For N = 3, ^[7] = 1 at 
the allowed pure-state density matrices defined by a dif¬ 
ferent circle 712 = (711 — 1.5)^ — 0.5^. Also at fV = 4 , 

F[7] = 2 at density matrix 711 = 2, 712 = 0. All these 
integer parts of the exact functional are pictured in the 
supplementary information. S6 


All Y come from a wavefunction 


Tr[Yl =4 



0 7,2 


0 ■l',2 


Figure S 6 : ^[ 7 ] for A = 0,1, 2,3,4 electrons 


Figure S5: 


H. Other ensembles for A = 1.5 electrons 


^GWA ^ [A(((iia((i2/3) + 7l((/)2a(/)i/3)] 

+g [clAifjjiafjjiP) -k (1 - cl)A{(j)2a(j)2f3)] 


In the consideration of fractional numbers of electrons 
the argument of convexity of A vs A is often used to sim¬ 
plify the ensembles that have to be taken. If convexity 
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Ensemble of N=1 and N=2 



Figure S7: Different ensemble formation of -^’[ 7 ^ com¬ 
bining pure state wavefunction for = 0,1, 2, 3 electrons 

is true, the lowest energy ensemble will always be given 
by the combination of the two integers at either side, 
e.g. Tn^s = (1 ~ (5)rAr + (5rAr+i. However, convexity 
has not been proven, with definitely known counterex¬ 
amples for certain electron-electron interactions, which 
indicates that most certainly convexity is not a general 
property of Hamiltonians [10]. Here, we test convexity 
for the two-site Hubbard hamiltonians, by taking ensem¬ 
bles of different electron numbers. We consider different 
pair-wise ensembles Fat^i.s = -|- 6|^'„2)(^'„2| 

with {ui,712} = {1,2}, {0,2}, {1,3}, {0,3}. We have also 
considered all possible ensembles, including those of three 
and four different particle numbers up to TV = 4, all of 
these lie higher in energy. 
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I. Supplementary Animations 

Supplementary animated gifs can be found in the 
arXiv source file or currently available via the following 
hyperlinks 

1) Varying t with fixed Ae 

2) Electron transfer by varying Ae for t = 1.0,0.2, 0.05 

3) Varying t for fractional number of electrons, N 







